function [J grad] = nnCostFunction(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda)
    Theta1 = reshape(nn_params(1:hidden_layer_size * (input_layer_size + 1)), hidden_layer_size, input_layer_size + 1);
    Theta2 = reshape(nn_params(hidden_layer_size * (input_layer_size + 1) + 1:end), num_labels, hidden_layer_size + 1);

    m = size(X, 1);
    E = eye(num_labels);

    for i = 1:num_labels
        Y0 = find(y == i);
        Y(Y0, :) = repmat(E(i, :), size(Y0, 1), 1); % Y 5000*10
    end

    X = [ones(m, 1), X];
    a2 = sigmoid(X * Theta1'); %5000*25
    a2 = [ones(m, 1), a2];
    a3 = sigmoid(a2 * Theta2'); %5000*10

    temp1 = [zeros(size(Theta1, 1), 1), Theta1(:, 2:end)];
    temp2 = [zeros(size(Theta2, 1), 1), Theta2(:, 2:end)];
    temp1 = temp1.^2;
    temp2 = temp2.^2;

    cost = Y .* log(a3) + (1 - Y) .* log(1 - a3); %cost 5000*10
    J = -1 / m * sum(cost(:)) + lambda / (2 * m) * (sum(temp1(:)) + sum(temp2(:)));

    delta_1 = zeros(size(Theta1));
    delta_2 = zeros(size(Theta2));

    for t = 1:m
        %step 1
        a_1 = X(t, :)';
        z_2 = Theta1 * a_1;
        a_2 = sigmoid(z_2);
        a_2 = [1; a_2];
        z_3 = Theta2 * a_2;
        a_3 = sigmoid(z_3);
        %step 2
        err_3 = zeros(num_labels, 1);

        for k = 1:num_labels
            err_3(k) = a_3(k) - (y(t) == k);
        end

        %step 3
        err_2 = Theta2' * err_3;
        err_2 = err_2(2:end) .* sigmoidGradient(z_2);

        %step 4
        delta_2 = delta_2 + err_3 * a_2';
        delta_1 = delta_1 + err_2 * a_1';
    end

    Theta1_temp = [zeros(size(Theta1, 1), 1), Theta1(:, 2:end)];
    Theta2_temp = [zeros(size(Theta2, 1), 1), Theta2(:, 2:end)];

    Theta1_grad = 1 / m * delta_1 + lambda / m * Theta1_temp;
    Theta2_grad = 1 / m * delta_2 + lambda / m * Theta2_temp;

    grad = [Theta1_grad(:); Theta2_grad(:)];
end
